function [Tn, cv] = as10_test_market_corr(m, N_dev, m_boot, g_ar, g_ar_boot, ...
    use_scope, brewer_m, N_brewer_dev, brewer_m_boot, brewer_g_ar, brewer_g_ar_boot,...
    est_corr, corr_m, N_dev2, corr_m_boot, corr_g_ar, corr_g_ar_boot, ...
    N_boot, reps, sig_level)
% purpose: compute the test statistic and the critical value

seed = RandStream('mt19937ar','Seed',100);
small_cons = 1e-10;

% N_obs: number of product-market obs; N_ineq: # of inequalities
[N_obs, N_ineq] = size(m); % Each N_dev rows correspond to a market.    
n_iv = size(g_ar, 2);
N_mkt = N_obs/N_dev;

%% compute the moments and the test statistic
miv = nan(N_mkt, N_ineq, n_iv);
miv_boot = nan(N_mkt, N_ineq, n_iv, N_boot);

for a = 1 : n_iv
    % interact the moments with instruments
    miv_a = m.*g_ar(:, a);%N_obs x N_ineq

    % average across potential products within a market: N_mkt x N_ineq
    miv(:, :, a) = mean(permute(reshape(miv_a', [N_ineq, N_dev, N_mkt]), [3 1 2]), 3); 

    if ~isempty(m_boot) %repeat for the bootstrap sample
        % interact the moments with instruments
        miv_boot_a = m_boot.*g_ar_boot(:, a, :); % (N_mkt*N_dev) x N_ineq x N_boot

        % average across potential products within a market: N_mkt x N_ineq x N_boot
        miv2_boot_a = permute(reshape(permute(miv_boot_a, [2 1 3]), [N_ineq N_dev N_mkt N_boot]), [3 1 4 2]);
        miv_boot(:, :, a, :) = mean(miv2_boot_a, 4);
    end
end

if use_scope 
    n_brewer_iv = size(brewer_g_ar, 2);
    brewer_miv = nan(N_mkt, N_ineq, n_brewer_iv);
    brewer_miv_boot = nan(N_mkt, N_ineq, n_brewer_iv, N_boot);

    for a = 1 : n_brewer_iv
        % interact the moments with instruments
        brewer_miv_a = brewer_m.*brewer_g_ar(:, a);
        % average across potential products within a market
        brewer_miv(:, :, a) = mean(permute(reshape(brewer_miv_a', [N_ineq, N_brewer_dev N_mkt]), [3 1 2]), 3);

        if ~isempty(brewer_m_boot)%repeat for the bootstrap sample
            % interact the moments with instruments
            brewer_miv_boot_a = brewer_m_boot.*brewer_g_ar_boot(:, a, :);

            % average across potential products within a market
            miv2_boot_a = permute(reshape(permute(brewer_miv_boot_a, [2 1 3]), [N_ineq N_brewer_dev N_mkt N_boot]), [3 1 4 2]);
            brewer_miv_boot(:, :, a, :) = mean(miv2_boot_a, 4);
        end
    end
else
    brewer_miv = [];
    brewer_miv_boot = [];
    n_brewer_iv = 0;
end

if est_corr
    n_corr_iv = size(corr_g_ar, 2);
    corr_miv = nan(N_mkt, N_ineq, n_corr_iv);
    corr_miv_boot = nan(N_mkt, N_ineq, n_corr_iv, N_boot);

    for a = 1 : n_corr_iv
        % interact the moments with instruments
        corr_miv_a = corr_m.*corr_g_ar(:, a);

        % average across potential products within a market
        corr_miv(:, :, a) = mean(permute(reshape(corr_miv_a', [N_ineq, N_dev2 N_mkt]), [3 1 2]), 3);

        if ~isempty(corr_m_boot)%repeat for the bootstrap sample
            % interact the moments with instruments
            corr_miv_boot_a = corr_m_boot.*corr_g_ar_boot(:, a, :);

            % average across potential products within a market
            miv2_boot_a = permute(reshape(permute(corr_miv_boot_a, [2 1 3]), [N_ineq N_dev2 N_mkt N_boot]), [3 1 4 2]);
            corr_miv_boot(:, :, a, :) = mean(miv2_boot_a, 4);
        end
    end
else
    corr_miv = [];
    corr_miv_boot = [];
    n_corr_iv = 0;
end

miv = cat(3, miv, brewer_miv);
miv = cat(3, miv, corr_miv);
M_mkt = miv(:, :); % N_mkt x (N_ineq*(n_iv+n_brewer_iv+n_corr_iv))

mean_M = mean(M_mkt);
mean_M_violate = mean_M.*(mean_M>0);

cov_M = cov(M_mkt);

Tn = sum((mean_M_violate.*1./(small_cons + sqrt(diag(cov_M)'))).^2);

%% critical value
if ~isempty(m_boot)% whether to compute the critical value as well 
    miv_boot = cat(3, miv_boot, brewer_miv_boot);
    miv_boot = cat(3, miv_boot, corr_miv_boot); % N_mkt x total_#_of_inequalities x n_iv x N_boot
    miv_boot = reshape(miv_boot, [N_mkt N_ineq*(n_iv+n_brewer_iv+n_corr_iv) N_boot]);

    % use the simulated var profits to estimate the sample variance of the moments
    M_mkt_boot = permute(miv_boot, [2 1 3]); %N_ineq*(n_iv+n_brewer_iv+n_corr_iv) x N_mkt x N_boot
    M_mkt_boot = M_mkt_boot(:, :)'; % (N_mkt*N_boot) x N_ineq*(n_iv+n_brewer_iv+n_corr_iv)

    Sigm = cov(M_mkt_boot) + small_cons;
    Dinvsqrt = 1./sqrt(diag(Sigm));
    Omg = diag(Dinvsqrt)*Sigm*diag(Dinvsqrt);
    Omg = Omg*1/2 + Omg'*1/2 + diag(diag(Omg))*(small_cons);

    % GMS selection
    mean_M = mean(M_mkt);

    % Andrews and Soares 2010 (4.5)
    xin = 1/sqrt(log(N_mkt))*sqrt(N_mkt)*Dinvsqrt'.*mean_M.*(mean_M<0); % 1 x N_moments
        
    sim_Z = randn(seed, reps, N_ineq*(n_iv+n_brewer_iv+n_corr_iv))*chol(Omg) + xin;
    Tn_S2_reps = sum(((sim_Z>0).*sim_Z./sqrt(diag(Omg)')).^2, 2);
    cv = quantile(Tn_S2_reps, 1 - sig_level)/N_mkt;
else
    cv = 0;
end







